Searching for topological density wave insulators in multi-orbital square lattice systems 
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We study topological properties of density wave states with broken translational symmetry in two-dimensional 
multi-orbital systems with a particular focus on t2g orbitals in square lattice. Due to distinct symmetry properties 
of d-orbitals, a nodal charge or spin density wave state with Dirac points protected by lattice symmetries can 
be achieved. When an additional order parameter with opposite reflection symmetry is introduced to a nodal 
density wave state, the system can be fully gapped leading to a band insulator. Among those, topological density 
wave (TDW) insulators can be realized, when an effective staggered on-site potential generates a gap to a pair 
of Dirac points connected by the inversion symmetry which have the same topological winding numbers. We 
also present a mean-field phase diagram for various density wave states, and discuss experimental implications 
of our results. 

PACS numbers: 



I. INTRODUCTION 

Identifying topological insulators has been one of the most 
fascinating research fields in contemporary condensed mat- 
ter physics. '"'' Topological insulators have a bulk gap like 
band insulators, but are distinguished by topologically pro- 
tected conducting edge states preserving time-reversal invari- 
ance. In particular, two-dimensional topological insulators are 
known as quantum spin Hall insulators with finite counter- 
propagating spin currents on the edge, analogous to quantum 
Hall states. Haldane^ proposed that the fictitious magnetic 
fluxes in the honeycomb lattice lead to the quantum anoma- 
lous Hall insulator (or Chern insulator). Generalizing Hal- 
dane's model including time reversal invariant spin-orbit cou- 
pling, it was theoretically shown that such a quantum spin Hall 
insulator can exist in graphenei^ A two-dimensional semi- 
conductor system with a uniform strain gradient was also pro- 
posed to be a candidate.^ Later, the predicted edge states in 
HgCeTe quantum well systems^ were experimentally verified 
which confirmed the existence of two-dimensional topologi- 
cal insulators.^ 

The topological insulators in these systems normally exist 
due to strong spin-orbit coupling.^''" When the spin-orbit cou- 
pling preserves spin rotational symmetry about an axis, the 
counter-propagating edge modes which carry opposite spin 
quantum numbers result in quantum spin Hall insulators. It 
was shown that these modes are protected by time reversal 
symmetry even in the absence of spin rotational invariance. '*' 
It was further pointed out that an effective spin-orbit coupling 
term can be generated by spontaneous spin rotational symme- 
try breaking in an extended Hubbard model on the honeycomb 
lattice Ji In these studies, the structure of the honeycomb lat- 
tice plays an important role, as the tight binding model on 
this lattice possesses two Dirac points at the Brillouin zone 
corners. Therefore in low energy description, various gapped 
insulating phases proximate to the Dirac semi-metal can be 
understood in terms of mass perturbations to gapless Dirac 
particles. For instance, the fictitious magnetic fluxes intro- 
duced by Haldane generate a mass term that has the opposite 
signs at the two Dirac points leading to an insulator with finite 



quantized Hall conductivity. The Dirac Hamiltonian approach 
further provides a framework to understand the time-reversal 
invariant Z2 topological insulators.'" 

While systems on the honeycomb lattice such as graphene 
naturally support two-dimensional massless Dirac particles 
in the bare band structures, this is not the case in a simple 
square lattice system which is an effective model for abun- 
dant layered perovskite materials in nature. In this respect, 
it is interesting to note that the recently proposed nodal den- 
sity wave state'^ exhibits gapless Dirac particles via broken 
translational symmetry. This proposal was made in the con- 
text of iron pnictide systems, where d-orbitals of t2g bands in 
an effectively two-dimensional square lattice give rise to sev- 
eral Fermi pockets with interesting topological properties. In 
this system, the spin density wave instability with the finite 
ordering wave vector Q = (tt, 0) (or (0, tt)) leads to band 
touchings between the k states with the momentum differ- 
ence of Q. In general, the degeneracies at the band touching 
points disappear because of the finite overlap matrix between 
the degenerate states induced by the density wave order pa- 
rameter However, in multi-orbital systems, because of the 
distinct symmetry properties of orbitals, the degeneracies at 
some band touching points are protected leading to nodal den- 
sity wave states, which is generally valid for any density wave 
orders. 

In this work, we ask if topological insulators can be 
emerged by gapping nodal points turning the system from 
nodal density wave states to topological density wave (TDW) 
insulators. To find such a TDW insulator, we first inves- 
tigate the properties of the nodal density wave states. We 
find that one general and important characteristic of the Dirac 
particles in nodal density wave states is that a pair of Dirac 
Hamiltonians connected by the inversion symmetry have the 
same topological winding numbers. Thus an effective stag- 
gered on-site potential generating a mass term, which has the 
same signs at the inversion symmetric nodal points, induces 
TDW insulators. This can be contrasted with the topologi- 
cal properties of the Dirac particles in the honeycomb lattice 
where the Dirac Hamiltonians at the two inversion symmetric 
nodal points have the opposite winding directions. '■^•''* Thus 



the mass term induced by, for example, a staggered sublattice 
chemical potential, which has the same signs at the two Dirac 
points would generate a topologically trivial band insulator as 
shown in graphene system. ^'^ 

The rest of the paper is organized as follows. In Sec. II, 
we first consider a simple two-band model Hamiltonian com- 
posed of dxz and dyz orbitals on the square lattice. After clas- 
sifying all possible charge and spin density wave order param- 
eters with the ordering wave vector Q — (vr, 0) based on their 
transformation properties under lattice symmetries, we estab- 
lish general relations between the locations of Dirac nodes 
and order parameter symmetries in Sec. III. The fact that d^z 
and dyz orbitals have the opposite eigenvalues under reflection 
symmetries along high symmetry directions in the momentum 
space, plays the key role for the emergence of Dirac points. In 
addition to the Dirac points coming from the Brillouin zone 
folding, additional contributions from quadratic band degen- 
eracy splitting are also discussed. In Sec. IV, topological 
properties of gapped density wave phases with two order pa- 
rameters with the opposite reflection symmetries are studied. 
Fully-gapped insulating phases can be obtained by introduc- 
ing two density wave order parameters which have the oppo- 
site eigenvalues under reflection symmetries. Among them, 
a certain combination turns the system to a TDW insulator. 
In Sec. V, the mean field phase diagram including the TDW 
phase is presented, which is obtained by solving an extended 
Hubbard model Hamiltonian with orbital degeneracy. Topo- 
logical density wave states in three orbital systems are dis- 
cussed in Sec. VI. Straightforward extension to three-orbital 
systems shows the general applicability of the idea we pursue 
in this work to obtain topological insulators in multi-orbital 
systems. Finally, we conclude in Sec. VII. 



n. TWO BAND HAMILTONIAN AND SYMMETRIES OF 
ORDER-PARAMETERS 

A. Tight-binding Hamiltonian 

We consider a tight binding Hamiltonian on the square lat- 
tice with two orbital (dxz, dyz) degrees of freedom at each 
site. A generic Hamiltonian which contains all possible hop- 
ping processes allowed by lattice symmetries is given by 



orbital states. In the above, 



^o = E4.^(k)^k,., 



k,CT 



(1) 



where 



i7(k) - (e+(k) - y.)l + e_(k)r3 + £.j,(k)ri. 



(2) 



Here a two-component field il\„ = [d^J.^ g-C^)' '^Iz o-C^)] de- 
scribes the creation of particles with d^z and dyz orbital fla- 
vors with spin a, and the Pauli matrix r connects these two 



e_i_(k) = — (ii + i2)(cosfc^ + cosky) — 4^3 cosfc^, cosfc, 
e_(k) = — (ii — t2){coskx — cosky) 
Exyiii-) = — 4i4 sinfca; sinfcj,. 



y 



yh 



(3) 

Diagonalization of i/(k) gives rise to the following two 
band dispersions. 



E±{\) = £+(k) - A* ± ./e2_(k) + 4 (k). 



(4) 



In addition to time-reversal symmetry T, the Hamiltonian 
iJo has the C4 point group symmetry, which consists of the 

four-fold rotation Cil, the inversion /, and the two reflec- 

2 

tions Px and Py mapping x to —x and y to —y, respectively. 
Each symmetry operation transforms a two-component field 
^a(kx,kij) in the following way. 
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Px 

Py 
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'ipaikx,ky) ^iT2'lpai-ky,kx), 



y^ayf^x^ f^y 



-^ 



'T3^<j{-kx,ky), 



Ipaikx^ky) -^ T3^pa{kx, -ky), 
4'tT{kx,ky) -^ -Ta'4)a{-kx,-ky). 



(5) 



If we choose the hopping parameters in such a way as ti=- 
1.0, ^2=1 -3, is = ^4 = -0.85, the i/o works as an effective 
two-band Hamiltonian describing the Fe-pnictide systems.— 
Given the hopping parameters above, the Fermi surface con- 
sists of two hole pockets and two electron pockets when the 
system is near half-filling.'-''^ A pair of electron and hole 
pockets are connected by a nesting wave vector Q — (tt, 0) 
(or (0, tt)), which drives various density wave instabilities.— 
Here we choose Q = (tt, Q)^ and perform a detailed study 
about the band structures of density wave ground states con- 
sidering all possible density wave order parameters. 



B. Symmetry of order parameters 

We consider various on-site density wave order parameters 
and investigate their symmetry properties. Since we have two 
orbitals per site, there are 4 different on-site charge density 
wave (CDW) states with the ordering wave vector Q = (tt, 0), 
which are given by 



1 ^ 

^* = ]^EE dlA^)[nUd,A^^ 



(6) 



k a, 6=1 



where a, h are indices describing the dxz (a=l) or dyz (a=2) 
orbital states. N counts the number of unit cells in the sys- 
tem. Similarly, we also define spin density wave (SDW) states 
choosing the spin ordering direction along the z-axis. 



(a)D3-CDW 





FIG. 1: (Color online) Description of representative density wave 
ordering patterns with the ordering wave vector Q = (7r,0). (a) 
Di charge density wave (Da-CDW) ordering, djjz and dyz or- 
bitals align alternatively along the x direction while local charge 
and spin densities are uniform, (b) M2 spin density wave (M2- 
SDW) ordering. Here orbital and spin orderings occur at the same 
time. A solid circle represents the d+ orbital defined as d+ = 
{dxz + idyz)/V2 while a dotted circle indicates the d_ orbital given 
by d_ — {dxz — idyz)/\^. The arrows inside circles describe spin 
ordering. 



k a,b=l 

These 8 order parameters represent distinct phases with dif- 
ferent broken symmetries. For example, the D^ CDW order 
parameter corresponds to^j(—l)'^(ni.a;z — ni^j,z) where ni,a 
is the density of electrons with the orbital a at the site i. Thus 
it is characterized by the relative density difference between 
two orbitals (orbital ordering), which alternates along the x di- 



rection, while keeping the total density {ni^^z 



''i.yz 



) uniform 



on every site as shown in Fig.[T](a). This breaks translational 
symmetry doubling the unit cell along x direction. On the 
other hand, the A/2 SDW order parameter described in Fig.[T] 
(b) corresponds to a staggered spin-orbit coupling. This is be- 
cause M2 can be written as X]i(~l)*^'5'i,zii,z where Li_z is 
proportional to [di^xz + idi,yz)Hdi^xz + idi,yz) - {di,xz - 



^^i,yz) y^i^xz 



H,yz 



2*^i ajz'^i.yz + ^■^■- However, un- 



like the uniform spin-orbital coupling J^i Si^zLi^z, the M2 is a 
staggered spin-orbit coupling with alternating signs along the 
X direction. It breaks spin-rotational and translational sym- 
metries but preserves time reversal symmetry. In addition, _Do 
and A/q describe conventional charge and spin density wave 
states, respectively. It was found that A/q describes the lead- 
ing density wave instability in Fe-pnictides.'- 

The above 8 order parameters can be distinguished by their 
transformation properties under lattice symmetries. The sym- 
metries of density wave order parameters are summarized in 
Table H] Note that every density wave state has even parity un- 
der the inversion symmetry. Moreover, all the diagonal den- 
sity wave states Di (or Mi) with i = or 3 are even under the 
two reflection symmetries while the other off-diagonal den- 
sity wave states with i = 1 or 2 are odd under the reflections. 
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TABLE I: Symmetry of density wave order parameters. Here '+' ('- 
') indicates 'even' ('odd') symmetry of the order parameters under 
the corresponding symmetry operation. 



These symmetry properties of density wave order parameters 
strongly constrain the location of Dirac nodes generated by 
the Brillouin zone folding and the winding numbers around 
Dirac nodes in the momentum space, which are discussed in 
detail in the following section. 



III. NODAL DENSITY WAVE PHASES 

One intriguing property of the Q = (tt, 0) density wave 
ground states is that a large number of Dirac nodes emerge in 
the band structure.'^ The numbers and locations of the nodal 
points depend on band dispersions and the symmetries of the 
order parameters. 

There are two different sources generating nodal points in 
general. One way is via introducing a density wave order pa- 
rameter carrying a finite momentum. This induces a Brillouin 
zone folding which generates several band touching points. 
In most cases, the degeneracy at the band touching point is 
lifted because the density wave order parameter induces a fi- 
nite overlap between the pair of states touching at a point. 
Henceforth a band gap opens up. However, when the band 
touching occurs at a high symmetry point in the Brillouin 
zone, the overlap matrix vanishes due to the lattice symme- 
tries, generating symmetry protected nodal points. 

The second group of nodal points come from the splitting of 
quadratic band touching points, which exist in the bare band- 
structure. Because of the underlying four-fold rotational sym- 
metry, the original hopping Hamiltonian in Eq.([ril supports 
quadratic band crossing points. '21^-22 -pj^g introduction of the 
density wave order parameter carrying a finite momentum 
splits a quadratic band touching point into two Dirac points 
along high symmetry directions in the momentum space. In 
the following, we discuss in detail the relation between the 
order parameter symmetry and the locations of Dirac points 
derived from these two different sources in separate subsec- 
tions. 



A. Dirac nodes generated by Brillouin zone folding 



We first focus on the generation of Dirac nodes along the 
fcy-axis. In Fig.|2a) (Fig.|2jb)), we plot the energy dispersion 
of the two bands given in Eq. ^ along the fc^r = {k^ — tt) 
direction. Since the e^y term in Eq. ([3]l, which describes the 
hybridization between d^z and dyz orbitals, vanishes along 
the kx — Q axis, the upper and lower bands in Fig. |2a) are 
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FIG. 2: (Color online) (a) Band dispersion along k^ = for the 
hopping Hamiltonian in Eq.lUll. Here we use a red solid (blue dotted) 
line to indicate a band which is odd (even) under the Px reflection 
symmetry, (b) Band dispersion along the kx=n direction, (c) Band 
structure along kx=0 after the unit cell doubling due to the Q=(7r,0) 
density wave ordering. Four bands in (a) and (b) meet and disperse 
together along kx=0 after the Brillouin zone folding. Note that the 
zone folding generates 8 band touching points. Black solid (dotted) 
circles indicate the band touching points between two bands having 
the same (opposite) P^ eigenvalues. 



just dyzand dxz bands, respectively. For k^ = 0, H{k) is 
invariant under the P^ reflection symmetry which transforms a 
momentum k^ to —k^- Therefore each band is an eigenstate of 
Px with eigenvalues of ±1. This is consistent with the fact that 
dxz (dyz) orbital is odd (even) under Px- Similar analysis can 
also be applied to the two bands dispersing along the kx = tt 
direction. Since H{k) has Px symmetry along kx ~ tt, the 
two bands also have definite Px eigenvalues. In Fig.|2] the Px 
even (odd) bands are represented by blue dotted (red solid) 
lines. 

Once we introduce a density wave order parameter with the 
ordering wave vector Q — (tt, 0), the unitcell doubles along 
the x-direction, which leads to the Brillouin zone folding in 
the momentum space. Thus within the reduced Brillouin zone, 
we have four bands dispersing along the ky-axis. Note that 
the zone folding generates 8 band touching points, which are 
indicated by circles in Fig.|2|c). Here the band touching point 
between two bands with the same (opposite) Px eigenvalues 
is encircled by a solid (dotted) circle. 

The degeneracy between two states, |V'i(k)) and |^2(k + 
Q)), touching at the momentum k after the Brillouin zone 
folding, is lifted when the matrix element of the density wave 
order parameter Di between these two states is finite, that is, 
('0i(k)|£)i|V'2(k + Q)) ^ 0. Therefore if the order parameter 
Di (or Mi) is Px even, the degeneracy is lifted when the two 
degenerate bands have the same Px eigenvalues. However, the 
nodal point remains gapless if the two degenerate bands have 
the opposite Px eigenvalues. 




FIG. 3: (Color online) The band structures of the Q = (tt, 0) density 
wave ground states along the fcy-axis. (a) For P^ even density wave 
states, (b) For Px odd density wave states. 



On the other hand, if the density wave order parame- 
ter Di (or Mi) is Px odd, the full Hamiltonian is not in- 
variant under Px anymore. However, even in this case 
(-01.31 (k)l-Di 1-02.32 C^ + Q)) = in the weak coupling limit, 
if si — S2 where s refers to Px eigenvalues. Namely, the 
matrix element of Di, which is odd under Px, vanishes when 
the two degenerate eigenstates have the same Px eigenvalues. 
To understand this point clearly, let us define the eigenvector 
|(/'n,s(k)) of the hopping Hamiltonian Hq with the even (s = 
+) or odd (s = -) Px eigenvalues. Here n is a band index. Now 
we turn on a small density wave order parameter Di which 
is odd under Px- Since Px eigenvalue is not a good quantum 
number, |0„,s(k)) can be contaminated by the states with the 
opposite Px eigenvalue IcjylJgCk+Q)) , leading to 



*i°l 



(k))^|0„,3(k)) 



.(k)) 



.(k+Q)), 



where |(/)3(k)) — X^n c„|(/'n,s(k)) is a linear combination of 
the states with the Px eigenvalue of s, while {(pgik+Q)) = 

Sn '^'nl'f'n s(k+Q)) IS a linear combination of the states with 
the opposite Px eigenvalue of s. Notice that |03(k)) and 
|03(kH-Q)) have a momentum difference given by the order- 
ing wave vector Q carried by the density wave order parameter 
Di. Because of the fact that the two components of the wave 
function with the opposite Px eigenvalues have the momen- 
tum difference given by Q, it is straight forward to show that 
(i/;i,3i (k)| A|V^2,S2 (k + Q)) = if si - S2. 

Therefore the nodal point remains gapless if the order pa- 
rameter Di is Px even (Px odd) while the two generated bands 
have the opposite (same) Px eigenvalues. It means that 4 
nodal points among the 8 band touching points remain gap- 
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FIG. 4: (Color online) (a) Band dispersion along ky — centered at 
k=(0,0). (b) Band dispersion along ky = centered at k=(7r,0). (c) 
Band structure along ky=0 after the unit cell doubling. Four bands in 
(a) and (b) meet and disperse together along ky=0 after the Brillouin 
zone folding. Here we use a red solid (blue dotted) line to indicate 
a band which is odd (even) under the Py symmetry. Black solid 
(dotted) circles indicate the band touching points between two bands 
having the same (opposite) Py eigenvalues. 



less independent of the condition that the order parameter is 
even or odd under P^ reflection symmetry. 

In Fig. |3] we plot the band structure of the density wave 
ground states along the fcj,-axis. Fig.|3](a) corresponds to the 
density wave orders Do, D^, Mo, M^, which are P^ even, 
while Fig. |3](b) describes the band structure for the other or- 
der parameters, Di, D2, Mi, A/2, which are odd under P^ 
symmetry. Notice that nodal points show opposite behavior 
for these two different classes of order parameters. Namely, 
when a nodal points remains gapless for one order parameter, 
it is gapped out for the other order parameter which has the 
opposite Px eigenvalue. 

We can extend the same analysis to understand nodal points 
lying along the fc^^-axis. In this case we have Py reflection 
symmetry mapping y to — y. However, compared to the previ- 
ous analysis for nodal points on the fcj^-axis, there is one im- 
portant difference in this case. Before the unitcell doubling, 
we have two bands dispersing along the fc^r-axis. The Bril- 
louin zone folding induces overlaps of these two bands with 
themselves. In Fig.H) we plot the dispersion of the two bands 
along the fc^^-axis centered at k = (0, 0) (Fig. lUa)) and at 
k = (tt, 0) (Fig. |5b)). The 4 bands after the zone folding 
displayed in Fig. Uc) can be obtained by superposing the 4 
bands in Fig.Ua) and (b). In Fig.Hfc), we plot the bandstruc- 
ture from k^ = — tt to kx = tt for convenience although the 
first Brillouin zone is from kx = — 7r/2 to kx = tt/2. Note 
that in Fig. |4jc) the location of solid and dotted circles are 
interchanged compared to those in Fig. |2c). Because of this 
difference, the location of Dirac nodes along the kx and ky 
axes also show the opposite behaviors. 



B. Dirac nodes from quadratic band crossing 

The band structure of the two-band hopping Hamiltonian 
Hf) in Eq. ^ supports two quadratic band crossing points at 
k = (0,0) and k = (tTjTt).'^'^^ Splitting of these quadratic 
band crossing points generates additional Dirac points, which 
contribute additional Chern numbers for various insulating 
phases. 

Expanding the Hamiltonian H{k) in Eq. (|2]l near k = 
(0, 0), we obtain the following low energy effective Hamil- 
tonian, 



H, 



eff 



d^k^Hk)H^,,d{k)i;ik), 



(8) 



in which 

i?quad(k) = a{kl + fc2)fo + (3kxkyTi + -f{kl - kl)T3, (9) 

where a = (ii +i2+4i3)/2, /3 = -4^4, and 7 = (ii-t2)/2. 
Nontrivial topological property of the quadratic band crossing 
point is reflected in the winding number N^, which is defined 
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l/dk.(V;t(k)|Vk|7/.(k)}, (10) 

where |V'(k)) is a Bloch wave function corresponding to one 
of the bands involved in the band touching and C is a closed 
loop in the momentum space encircling the band crossing 
point. A quadratic band crossing point contributes A'^vK = i2, 
which is twice larger than the winding number around a Dirac 
point. '"riLSl 

Adding a generic perturbation given by F = J2i=i ^^i^i^ 
the degeneracy at the quadratic band crossing point can be 
lifted. ?7i2 term breaks time-reversal symmetry and the de- 
generacy is lifted by opening a gap. On the other hand, rrii 
and TO3 terms that break 4-fold rotational symmetry, split the 
quadratic band touching point into two Dirac pointst^^ 

Now we consider the effect of the Q = (tt, 0) density wave 
orderings on the degeneracy lifting at quadratic band cross- 
ing points. Since the density wave order parameters carry the 
momentum Q = (tt, 0), they cannot couple to the degenerate 
states at k = (0, 0) (or k — (tt, tt)) at first order. The low- 
est order contribution to degeneracy lifting at quadratic band 
touching points starts from the second order processes. We 
first consider charge density wave order parameters given by. 



HcDw = ^ V'k,^-D?/'k+Q,<T, 



(11) 



k,o- 



where D = X),;=o ^i^i- Treating the above Hcdw as a per- 
turbation, the standard second order perturbation theory gives 
rise to the following effective Hamiltonian near the quadratic 
band touching point at k = (0, 0), 



3 
=i?quad(k) + ^ mp'f,; 



(12) 



1=0 



(a) 



nodal points coming from the Brillouin zone folding induced 
by the Q — (tt, 0) density wave states. On the other hand, 
red dots result from the splitting of quadratic band touching 
points. Two quadratic band touching points generate four 
Dirac points lying on the y-axis. With the understanding of 
the origin and locations of Dirac points, below we discuss how 
to achieve TDW insulators. 
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FIG. 5: (Color online) Distribution of Dirac points for the Da (or 
M3) density wave state, (a) Four bands within the reduced Brillouin 
zone. The energy eigenvalue Ei (k) for a band i satisfies Ei (k) > 
E2{k) > Ealk) > E4{k). (b) Dirac points between the band 1 and 
2, (c) between band 2 and 3, (d) between band 3 and 4. Blue (Red) 
dots indicates Dirac points coming from the Brillouin zone folding 
(splitting quadratic band touching points). 



in which 



r4°) =X{{Dl + Dl){h +t2+ 4i3) 

+ (Z?o + D3f{t2 + 2*3) + [Do - D:if{h + 2h)}, 
m[,'' ^2\Di{D:i(-h + ^2) + D^ih + ^2 + 4^3)}, 
m^"' ^2\D2{D:i{-ti + tj) + Do{ti + ta + 4^3)}, 
m^ ^\{{Dl+Dl){h-t2) 

+ {Do + D3f{t2 + 2U) - [Do - D:if{h + 2h)}, 



where A == -l/{8(ti + 2t^){t2 + 2*3)}. Note that as long as 
only one of the order parameters has finite magnitude while 



all the other order parameters are zero, mp 



(1) 



,(2) 



0. 



In other words, if _D„ / for a given n while all the other 
Di-in — 0, the quadratic band crossing point always splits 
into two Dirac points along the main axes. 

Combining the contributions both from the Brillouin zone 
folding and from the splitting of quadratic band crossing 
points, we show the distribution of Dirac points for a Z?3 (or 
A/3) density wave ground state in Fig. |5] There are four bands 
within the reduced Brillouin zone as shown in Fig. EJa). We 
assign a band index i such that the energy eigenvalue Ei (k) 
of the band i satisfies ^i(k) > £'2(k) > ^3(k) > ^4(k). 
The location of Dirac points between the upper two bands 
(band 1 and 2) are indicated in Fig.|5jb). Similarly, the Dirac 
points between the middle (bottom) two bands are described 
in Fig. |5jc) (Fig. |5jd)). Notice that there are many Dirac 
touching points between the bands. Blue dots indicate the 



IV. TOPOLOGICAL PROPERTIES OF THE GAPPED 
DENSITY WAVE PHASES 



A single density wave order parameter induces a metallic 
phase with many Dirac points. The locations of Dirac points 
are determined by the transformation properties of the order 
parameters under the reflection symmetries P^ and Py. There- 
fore to get an insulating phase, two coexisting density wave 
states, in which one is even and the other is odd under the 
Px and Py symmetries, are required. In addition, according 
to the order parameter symmetries summarized in Table |I] if 
time-reversal invariance is imposed, there are only four differ- 
ent ways of choosing a pair of density wave order parameters, 
which give rise to a gapped phase. The four pairs of time 
reversal invariant density wave order parameters with the op- 
posite transformation properties under the reflections P^ and 
Py, are given by (i?3, Di\ (D3, A/2), (Do, /?i), and (Do, 
A/2). 

Since the z component of the spin, Sz is conserved, the 
Chern number N^'^ is well defined for each band in a fully 

gapped phase. ^^"^^ Here N^l (N^!) is the Chern number of 
the nth spin-up (spin-down) band. For every pair of the den- 
sity wave order parameters generating a fully gapped phase, 
the four bands within the reduced Brillouin zone are well- 
separated from each other with a finite gap between any pairs 
of the bands. Each band is distinguished by the index n rang- 
ing from 1 to 4 as the energy decreases. The Chern number 

r(") 



(13) ^c (T °f ^^^ ^^^ band with the spin a is defined as. 



N, 



(n) 
C,cr 



= ^ [ d'kF^^'\k), (14) 
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where the momentum space Berry 
for the 71th band with the spin 



Fi")(k) E 

the Berry potential j4/j"i^(k) 



dk^4:Uk) 



curvature Fa (k) 
a is defined as 

in which 



given by A}i"^(k) 



~i{^r' {k)\dk^^ |$J,"^(k)),2L29 Here the Bloch wave function 
$CT (k)) is defined within the reduced Brilluin zone (RBZ). 
Explicit computation of the Chern numbers using Eq. (fl4l i 
shows that every band of the insulating density wave phase 
with finite D3 and A/2, has a nonzero Chern number as shown 
in Table |II] On the other hand, every band has zero Chern 
number for the other three gapped phases defined with a pair 
of nonzero order parameters given by (£'3, Di), {Dq, Di), 
and (Dq, A/2). 





Net 


Nc,i 


band 1 
band 2 
band 3 
band 4 


+ 1 
-3 
+ 3 
- 1 


- 1 
+ 3 
-3 
+ 1 



TABLE II: The Chem number of each band for the gapped density 
wave ground state with finite D^ and A'h . 
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(16) 



Here the momentum q of the Hamiltonian ^^^'^'^(q) is mea- 
sured with respect to the degeneracy point at k = k^ (/i = 
+, — ). The two-component fermion field ^p.£r(q) is given by 



(a) Dsonly 

E(k) 



-n 



(b) D3+M2 
E(k) 





FIG. 6: (Color online) (a) Band structure of the D3 nodal density 
wave state along the ky axis, (b) Band structure of the gapped density 
wave phase with finite D3 and M2. 



*p.^(q) = 



aidy^^^{kf, + q) + a2d.y^,CT(kp + Q + q) 
I3id^z,a{^f, + q) + I32dxz,a{^f, + Q + q) 



(17) 
where the constant coefficients Ui and /3i satisfy of + a^ — 
/?^ + /3| = 1. Explicitly, a^ and ^i {i=l,2) are given by 



ai = (2 + ^A + DD/Js + 2Dl + 4^4 + 7^1, 



"2 = D^l 



/3i = (2 - 



P2 = ^3/ 
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Topological properties of the topological density wave 
ground state with finite D3 and M2 



Nontrivial topological properties of the fully gapped den- 
sity wave phase with nonzero D^ and M2 can be understood 
in the following way. We first consider the charge density 
wave state with the finite D^ order parameter The D^ density 
wave phase supports many Dirac points whose distribution is 
described in Fig.|5] Now we turn on a small A/2 which induces 
gap opening at each nodal point, leading to the fully gapped 
insulating phase, which is described in Fig.|6] The degeneracy 
lifting at each nodal point can be understood as a result of the 
mass perturbation, induced by the finite A/2, to the gapless 
Dirac particles. The nonzero Chern number of each band is 
obtained by adding up the Chern number contributions of the 
massive Dirac particles derived from the corresponding band. 

We first focus on the two Dirac nodes lying along the ky 
axis between the band 3 and 4 shown in Fig. EJd). The ef- 
fective Dirac Hamiltonian can be obtained by linearizing the 
Hamiltonian near the two nodal points sitting at the momen- 
tum k=k_|_ and k_. To simplify the computational procedures, 
the hopping parameters ti are slightly shifted from the ini- 
tial values given in Sec. Ill Al to ti — —t2 = H = ti=-\.Q. 
This small parameter change does not affect the topological 
properties of the gapped phase but shifts the nodal points to 
k± = (0, ±^) making analytical analysis simpler. The ef- 
fective Hamiltonian describing the low energy fermions near 
these two nodal points is given by 



H Dirac 
eff 



E 



/*l.,.(qX'r(q)*M.^(q) 



(15) 



Notice that the first (second) component of '^ ^_a (q) is derived 
entirely from the dyz (d^z) orbital. 

Now we include the A/2 spin density wave order parameter 
which generates a mass term in the low energy limit given by 



Hmass 
Ah 
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d'lHA^in- 



2\D.\M2a 



/?3v/4" 
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-^3 



fa}* 



p.<T(q)- 

(18) 



For the given spin a, the mass term has the same magni- 
tude and sign at the two Dirac points. At each Dirac point, 
this mass term opens a gap and contributes to the Chern 
number Nc.a — + 2\m\ for the upper band (band 3) and 
Nc.„ = -^^ for the lower band (band 4).24.30-32 Adding 



1C,a 



2 1 Mo 



the Chern number contributions from the two Dirac points, 
the total Chern number of the band 4 with the spin a is given 
by N^ [^ — — sgn(A/2)(T, which is consistent with the result 
obtained from the integration of the BeiTy curvature over the 
reduced Brillouin zone using Eq. ( fT4l l. (See Table [ill ) In the 
case of the band 3, the Chern number is determined after in- 
cluding the additional contributions from the nodal points be- 
tween the band 2 and band 3. 

Similarly, the trivial topological property of the gapped 
phase with finite D^ and /?i can also be understood by apply- 
ing the same analysis. For the D^ nodal density wave state, 
the small Di charge density wave order parameter generates 
mass perturbations to Dirac particles, which can be described 
by the following Hamiltonian, 



Tjmass 
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d\^lA<\){ 



2\D.\D, 
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n}*,..(q). 



(19) 
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Band filling 


Spin Chem Number (Cs) 


^2 invariant (u) 


! 
? 

4 


-2 
+ 4 
-2 


1 

1 



TABLE III: Spin Chem numbers and topological Z2 invariants for 
the gapped density wave ground state with finite D3 and M2. 



Notice that this term just induces the shifting of the nodal 
points away from the ky-axis. Once a Dirac point moves away 
from the reflection symmetry axis, the degeneracy at the band 
touching point is no longer protected by the symmetry and a 
gap opens, because the density wave order parameters support 
finite matrix elements between the two bands touching at the 
nodal point. Since the Di charge density wave order param- 
eter does not generate a mass term to the Dirac Hamiltonian, 
it has no contribution to the Chern number leading to the zero 
Chern numbers of all bands. 

We apply similar analysis to every Dirac point derived from 
the Brillouin zone folding for all pairs of density wave order 
parameters generating fully gapped phases. In all cases, it is 
confirmed that the Chern number of each band obtained by 
summing up the Chern number contributions from the Dirac 
points is identical to the result obtained by the integration of 
the Berry curvature in the momentum space. 

In addition, the coexisting D3 and M2 density wave order 
parameters also lift the degeneracies of the two quadratic band 
touching points at k = (0, 0) and (0, vr) leading to a fully 
gapped bandstructure. In contrast to the case of Dirac points, 
the Chern number obtained by lifting a quadratic band degen- 
eracy is two times larger than the contribution from a single 
Dirac point. However, since the two quadratic band crossing 
points lead to the Chern number contributions with the oppo- 
site signs, the net effect of the two quadratic band touching 
points vanishes. 

Since the four bands are well separated from each other for 
the topological density wave phase with D3 ^ Q and A/2 ^ 0, 
if the magnitude of the order parameter AI2 is large enough, an 
insulating ground state is obtained whenever the Fermi level 
lies in the gap between two neighboring bands. Therefore 
there are three different insulating phases, in principle, when- 
ever the Fermi level lies between the band n and n + 1 (n= 1 , 2, 
3) corresponding to the band filling factor iVfiiiing — (4 — n)/4. 
The topological property of the insulating phase can be explic- 
itly characterized by computing topological invariants. We 
first consider the spin Chern number Cs which is defined in 
the following way. 



Cs^ 
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(20) 
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where the summation includes all the occupied bands. When 
the Sz is conserved, the spin Chern number Cs is quan- 
tized and characterizes the two dimensional topological insu- 
lators.^^ In Table. Hn] we show the spin Chern numbers for the 
insulating phases. It is interesting that the spin Chern numbers 
are nonzero for all cases. Therefore as long as the z compo- 
nent of the spin is conserved, we can obtain the topological 



insulator with finite spin hall conductivity for every quarter 
filling. However, in the presence of spin non-conserving per- 
turbations and disorders, the spin Chern number is not well- 
defined and conserved only modulo 4?"^'^^ In other words, 
when the Sz is not conserved, the half-filled system is equiva- 
lent to the phase with zero spin Chern number, which is noth- 
ing but a topologically trivial band insulator However, it is 
important to notice that even in this case the system remains 
as a topological insulator with nonzero Z2 topological invari- 
ant for 1/4 and 3/4 fillings. 

We also compute the Z2 topological invariant v shown in 
Table |III1 which can be used to distinguish topological in- 
sulators (ly — 1) from trivial band insulators (ly — 0) for 
generic time-reversal invariant systems. Since the system has 
the inversion symmetry, the Z2 invariant can be obtained from 
the parity eigenvalues ^mi^i) of the occupied bands with 
the band indices m at the time-reversal invariant momenta 
Fi.^"^ Using the reciprocal lattice vectors G^ {i=l, 2), the four 
time reversal invariant momenta can be written as ri=nj^n2 
=(niGi + n2G2)/2 with ni,2 = 0, 1. Explicitly the Z2 topo- 
logical invariant ly is given by 



rii— 0,1 m 



(21) 



where the parity eigenvalues of the occupied bands at the four 
time reversal invariant momenta are multiplied. Because of 
the time reversal symmetry, each band is doubly degenerate 
at the time reversal invariant momentum and every Kramers 
doublet share the same inversion parity. The Z2 topologi- 
cal invariant counts the parity of one state for each Kramers 

■ 34 
pair 

As shown in Table, [nil the Z2 topological insulators exist 
when the band filling is one-quarter or three-quarter. How- 
ever, the 3/4 filled case requires unreasonably large AI2 to 
achieve an insulating phase. This is because, as shown in 
Fig. |6] the overall structures of the band 1 and 2 are in par- 
allel. To open a full gap between the upper two bands (band 1 
and 2), the magnitude of the M2 density wave order parameter 
should be as large as their bandwidth. Therefore the quarter- 
filled system is the most favorable for the realization of the Z2 
topological insulator 

To support further the nontrivial topological properties of 
the topological density wave phase with D3 ^ and M2 7^ 
0, we compute the edge state spectrum by considering the 
Hamiltonian on a strip geometry, which is infinite in the x- 
direction but finite in the y-direction with open boundaries at 
y=l and y=Ny. Here Ny indicates the number of lattice sites 
in the y-direction. The energy spectrum of the system with 
Ny = 40 is described in Fig.|2l which shows the existence of 
robust gapless edge states traversing between the lower two 
bands (band 3 and 4) and the middle two bands (band 2 and 
3). The upper two bands (band 1 and 2) are not well sep- 
arated for I?3 = 0.85 and M2 = 0.75, which are used to 
obtain the energy spectrum. Each edge state represented by a 
red dotted line is doubly degenerate, one with spin-up and the 
other with spin-down. For the 1/4-filling with the chemical 
potential lying in the gap between the band 3 and 4, there are 
two gapless edges states on each boundary propagating in the 




FIG. 7: (Color online) Energy spectrum for the topological density 
wave ground state with D3 — 0.85 and _A/2 = 0.75 in a strip ge- 
ometry. Here we use open boundary conditions with Ny =40 sites 
along the y-direction. The red dotted lines stand for edge states, all 
of which are doubly degenerate. For 1/2-filling, the number of edge 
states is twice larger than that for 1/4-filling. 




FIG. 8: (Color online) The winding directions of the d(q) vector 
around the Dirac point. Black arrows describe the two component 
d — {dx, dy) vector in Eq.l|23l) along the circular path around the 
Dirac point at the center, (a) For the two Dirac points between the 
band 3 and band 4 in the nodal density wave state with the finite Dz- 
The d vector has the same winding direction around the two Dirac 
points, (b) For the two Dirac points on the graphene system. The 
d vector has the opposite winding directions around the two Dirac 
points. 



opposite directions with the opposite spin quantum numbers. 
On the other hand, for the 1/2-filling, there are four gapless 
edge states on each boundary consistent with the fact that the 
spin Chern number Cs — 4 whose magnitude is twice larger 
than that for the 1/4-filling with Cs = —2. Therefore for 
the collinear spin ordering with M2 7^ which conserves Sz, 
there are robust gapless edge states for both the quarter-filled 
and half- filled systems. 



particles, respectively. 

The topological property of the Dirac particles in the D3 
nodal density wave state can be understood in the following 
way. The low energy Hamiltonian i/^'^'^ in Eq.dTSb for the 
Dirac particles near the momentum k = k±, where the band 
touching points between the band 3 and 4 locate, can be writ- 
ten as 



lff^'^(q)=/i±,,(q)fi+/i±,„(q)f3. 



(22) 



B. Comparison to the honeycomb lattice 

It is interesting that a simple on-site density wave order pa- 
rameter can generate insulating phases with nontrivial topo- 
logical properties. This result can be contrasted with the topo- 
logical insulator on the honeycomb lattice where complex sec- 
ond neighbor hopping processes are required to obtain a topo- 
logical insulator while the simple on-site staggered chemical 
potential gives rise to a trivial band insulator.'*'^'*' It is the dis- 
tinct topological properties of the Dirac particles in the D^ 
nodal density wave phase in the d-orbital system, which make 
it possible to realize the topological insulator by introducing a 
simple on-site order parameter M2. 

In this subsection, we discuss the topological property of 
the Dirac particles in the Z?3 nodal density wave state in de- 
tail and compare it with the topological property of the Dirac 
particles on the honeycomb lattice. In the forthcoming discus- 
sion we neglect the spin degrees of freedom and focus on the 
condition under which the insulating phase possesses a finite 
Chern number, which is nothing but a Chern insulator. Once 
we find the condition to obtain a Chern insulator, the time 
reversal invariant topological insulator can be realized by su- 
perposing two Chern insulators with spin-up and spin-down 



Since \ h\^{(\) + /ij_ (q) is nonzero away from the de- 
generacy point, a two component vector d± (q) with the unit 
length can be defined as 



d±{(\) = (d±,x(q),d±,y(q)) 



(/i±,x(q)>±,a(q)) 



hl^M + hl^yi^) 



(23) 



Along the circle Cr^ satisfying q^ + q^ = i?^ 7^ with the de- 
generacy point at the center, the 2D unit vector d± (q) defines 
a map from the circle Cr. to the unit circle S^. Since the fun- 
damental group TTi (5^) = Z, the 2D unit vector d± (q) has an 
integer-valued topological invariant, which is nothing but the 
winding number Ns,, defined in Eq.(fTOli. In terms of the 2D 
unit vector (i±(q), the winding number Ns^ can be rewritten 
as. 



1 r /-//7 

^w = 7-$ dez-{dx—) 
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where the loop integral is defined along the circle where the 
momentum q = Re^^^ 
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In Fig. [8ta) we describe the directions of the two- 
component d vector along the circular path around the Dirac 
points for the D^ nodal density wave state. Notice that the d 
vector has the same winding direction with the winding num- 
ber of Nvj=l around the two Dirac points. The two Dirac 
points share the same winding number because of the con- 
straint imposed by lattice symmetries. Since the first (second) 
component of the two-component Dirac field 'i'^{q) is given 
by the dyz (dxz) orbital state, ^^^(q) transforms to -\I'p(— q) 
under the inversion symmetry due to the odd parity of d^z and 
dyz orbitals. Here fl has the opposite sign of /i. This imposes 
the following constraint on the pair of Dirac Hamiltonians re- 
lated by the inversion symmetry, 



Hh 



'-{q) =i7°"-(-q). 



(25) 



This constraint guarantees the same winding numbers for the 
two Dirac Hamiltonians iji^uac ^j^^j ^pn-ac jj. j^ important to 
notice that every pair of Dirac Hamiltonians related by the in- 
version symmetry satisfies the same constraint for nodal den- 
sity wave phases. 

The fact that a pair of Dirac Hamiltonians related by the 
inversion symmetry have the same winding numbers is the 
distinct topological property of the Dirac particles in the 
nodal density wave states, distinguishable from the topolog- 
ical properties of the Dirac particles in the honeycomb lattice. 
In this system, the two Dirac points at the corners of the first 
Brillouin zone have the opposite winding directions, which is 
described in Fig. ISb). Since the inversion symmetry inter- 
changes the two sublattices of the honeycomb lattice, each of 
which comprises one component of the Dirac fermion field 
^^(q), the Dirac fermion field ^^^(q) transforms, for exam- 
ple, to Tx'i'fi{—q) under the inversion symmetry. This im- 
poses the following constraint to the two Dirac Hamiltonians 
related by the inversion symmetry. 



ff+""(q) -r,i/°'"=(-q)7 



(26) 



The additional Pauli matrix reverses the winding direction of 
one of the d vectors, leading to the two Dirac Hamiltoni- 
ans with the opposite winding numbers. Therefore these two 
Dirac points can be pair-annihilated when they are brought 
together by perturbations. '''■'"^ 

The relative winding numbers of the pair of the Dirac 
Hamiltonians related by the inversion symmetry, strongly con- 
strain the topological properties of the insulating phases ob- 
tained by mass perturbations to the Dirac particles. The in- 
troduction of a constant mass term H'^^^^ = mf2 to the Dirac 
Hamiltonian in Eq. (l22l i gives rises to the third component 
dz{q) of the corresponding d vector. '^'^^ Explicitly, for the 
massive Dirac Hamiltonian given by 



i/°''^'=(q) =/i,(q)fi + hy{q)h + mfs, 
the 3D unit vector d^^ is defined as 

(feD(q) ={dx{q),dy{q),dz{q)) 
_ {hx{q),hy{q),m) 



(27) 



hl{q) + hl{q)- 



(28) 



If we introduce, for instance, a positive mass term to the 
two Dirac Hamiltonians corresponding to the D^ nodal den- 
sity wave phase described in Fig.UJa), both of the ^30 vectors, 
which have positive z components, move along the northern 
hemisphere as the momentum q sweeps over the two dimen- 
sional momentum space. The net solid angles subtended by 
these two dan vectors over the entire momentum space are 
the same, each of which covers 2-k. At this point, it is useful 
to take into account the following relation between the Chern 
number of the valence band and the 3D unit vector d^^ for the 
two band Hamiltonian in Eq. (|27]l,'- 



Nr 



An 
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{dk^diD X dkydao)- 



(29) 



The above identity implies that the Chern number counts the 
number of times the 3D unit vector ^30 winding around the 
unit sphere over the Brillouin zone torus. Therefore when a 
constant mass term is added to the two Dirac points connected 
by the inversion symmetry, the Chern number Nq = ±1 if 
the two Dirac points have the same winding numbers, which 
is realized in the nodal density wave ground state. 

In contrast, the net solid angles covered by the two ^30 
vectors in the honeycomb lattice have the same magnitudes 
but with the opposite signs. Therefore the total solid angle 
covered by these two ^30 vectors vanishes. The vanishing 
Chern number contributions from the two Dirac points leads 
to the topologically trivial insulating phase when the simple 
mass term in Eq. dZTl i is introduced. This contrasting behav- 
ior of the pair of 3D unit vectors ^30 in these two systems 
results in the distinct topological properties of the insulating 
phases, the topological insulator in the d-orbital system and 
the topologically trivial band insulator in the honeycomb lat- 
tice when constant mass terms with the same signs are added 
to the pair of Dirac points connected by the inversion symme- 
try. However, in graphene if we introduce mass terms with the 
opposite signs at the two Dirac points, a topological insulator 
with a finite Chern number can be obtained, which is realized 
by considering the complex second nearest neighbor hopping 
processes on the honeycomb lattice. 



V. MEAN FIELD THEORY 

Now we address the question whether the TDW insulator 
with finite D3 and A/2 can be achieved in real systems. In par- 
ticular, by taking into account the interactions between elec- 
trons, we investigate the conditions to realize the TDW insula- 
tor via spontaneous symmetry breaking. Previous mean field 
studies on the two-orbital Hubbard model with the hopping 
Hamiltonian Hq in Eq. ([T]) show that the leading instability 
of the system is the uniform spin density wave phase (Mq- 
SDW) which is described by nonzero A/q. However, there ex- 
ist another density wave order parameters including imaginary 
charge and spin density wave states, which are competing with 
the uniform spin density wave state (A/q-SDW) with small en- 
ergy differences. '^'^^'■'^ Therefore if we include longer range 
interactions which are not included in the multi-orbital on-site 
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FIG. 9: (Color online) Mean field phase diagram. Here we set 
(7 = 5, Jh = and plot the ground state phase diagram varying 
the inter-orbital on-site repulsion U2 and the effective next nearest 
neighbor repulsion V2,rff = V2B — 2V2A- There is a finite range 
between D3 charge density wave phase (Da -CDW) and M2 spin 
density wave phase (M2-SDW) where the topological density wave 
phase (TDW) becomes the ground state. Solid (dotted) lines indicate 
the first (second) order phase transitions. 



Hubbard Hamiltonian, another competing ground states, for 
example, the topological density wave state (TDW), can be- 
come the ground state replacing the A/q-SDW state. 

Including on-site and inter-site electron-electron interac- 
tions, the full Hamiltonian is given by 

HfuW = Ho + ifonsite + -^intersite (30) 

in which 

i a^l,2 i 

i cri,cr2 

+ Jh Y^{4^^4^,id^,2Ad^a^ + H.c), (31) 

i 

and 



^intersite ^ViA ^ ^ ni^anj,a + ViB ^ '^i,l"'J:2 
(ij) a=l,2 (y) 

+ V2A ^ ^ n^.anj,a+V2B ^ "i,l"j,2, 
(fe» °=1.2 (fe» 

(32) 

where Hq indicates the hopping Hamiltonian in Eq.([T]i. For 
the on-site interactions described by i/onsite, the intra-orbital 
repulsion U, the inter-orbital repulsion U2, and the Hund's 
coupling Jh are considered. In the i^intersite describing 
the inter-site Coulomb interactions. Via (Vib) indicates the 
nearest-neighbor Coulomb repulsion between electrons in the 



same (different) kinds of orbitals. Finally, V2A (V2B) indi- 
cates the next nearest-neighbor Coulomb repulsion between 
electrons in the same (different) kinds of orbitals. 

To investigate the existence of the topological density wave 
state (TDW) with D3 ^ and A/2 / and its competition 
with the uniform spin density wave phase (Mq-SDW), we ap- 
ply a mean field approximation to the Hamiltonian i/fuii- The 
resulting mean-field Hamiltonian is given by 
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(33) 



in which 



eo^^{M^-Dl) 



^{M'2+^Dl) 



^iM^-Dl-Mi) + {V2A-^^)Dl (34) 



and 



- i2V2A - V2b)D3 + ^i-Ds + aMo), 



A, 



U 



{D3 + a Mo) 



U2 



D3 
Jh 



+ i2V2A~V2B)D3 + ^{D3 

Ai2,a- =^21,0- = —iaAh —iaM2. 



a Mo), 



(35) 



Note that the nearest neighbor Coulomb repulsions Via and 
ViB do not contribute to the mean field Hamiltonian be- 
cause the order parameters have the ordering wave vector 
Q=(vr,0). 

The order parameters Mo, M2, D3 are determined by solv- 
ing the following self-consistent equations. 



^r,2,o-"r,2 



dr,2,a), 



^0 = ^ E E(-i)"^<<i.^^--.i.- 

A^2 = -TT Yl E ("l)''"^(*4,l,<Trfr,2,<T - idl 2 Ar,i,<y)^ 



D. 



N 



^(r^.ry}cr=± 



Y Y.^-i)'-HdU,j^,i,.-di 



2,a-d^,2,a)- 



= (l'x 



)a=± 



(36) 



The chemical potential /i is also determined self-consistently 
to satisfy the half-filling condition. 

The resulting mean field phase diagram is shown in Fig. |9] 
Here we choose U — 5, Jh — and compute the ground 
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state phase diagram as a function of the inter-orbital on-site 
repulsion U2 and effective next nearest neighbor repulsion 
V2,eff = V2B — 2V2A- In the absence of the inter-site inter- 
actions V2.etf = 0, the uniform spin density wave phase (Mq- 
SDW) dominates the phase diagram consistent with the pre- 
vious studies. However, the inter-orbital on-site repulsion U2 
suppresses the uniform spin density wave states (A/q-SDW) 
which is diagonal in the orbital space, but promotes the M2- 
SDW, which is off-diagonal in the orbital space, to the ground 
state. On the other hand, the D^ charge density wave phase 
(D^-CDW) is strongly affected by the next nearest neighbor 
Coulomb repulsions Vza and ¥23- In particular, the V2A, 
the next nearest neighbor repulsion between the electrons in 
the same kinds of orbitals, strongly favors the I?3-CDW be- 
cause the staggered orbital ordering described by Ds-CDW 
can avoid the energy cost coming from V2a- Notice that there 
is a finite range in the parameter space where both M2-SDW 
and D3-CDW are nonzero realizing the topological density 
wave phase (TDW). 

The above mean field phase diagram is obtained for the half 
filled case where the topological property of the TDW phase 
is not robust against perturbations breaking Sz symmetry. On 
the other hand, the TDW insulator at 1/4 filling maintains its 
topological properties as long as time reversal symmetry is 
preserved. It was shown, in the study of the single orbital 
extended Hubbard model, that the next nearest neighbor inter- 
action {V2A for our model) stabilizes a stripe pattern charge 
ordering with the momentum Q = (tt, 0).^** This occurs at 
1/4 filling when the on-site Hubbard interaction U is much 
stronger than the hopping amplitude t satisfying U >> t, so 
that double occupancy is almost frozen. In our model, there 
are two orbitals of dxz and dyz- Similar to the single orbital 
case, we find that the next nearest neighbor interaction stabi- 
lizes Z?3 order Therefore, when the on-site intra- and inter- 
orbital interactions satisfy U, U2 3> t, the Z?3-CDW ordering 
should be favored at 1/4 filling. 

To get a TDW insulator, finite M2 is required in addi- 
tion to D^. As shown in Sec. II B, the M2 order order 
parameter is equivalent to a staggered spin-orbit coupling, 
^ J^ii^^Y'' '^i,z-^i,z- C>ne can show that when spin-orbit in- 
teraction is present, A/2 term can be induced as long as D3 
sets in, since D3 leads to unequal density between dxz and 
dyz orbitals. Therefore, we expect that the TDW insulator can 
be obtained by tuning V2B when U,U2 ^ t at 1/4 filling in 
the presence of the spin-orbit coupling. 



VI. TOPOLOGICAL INSULATORS IN THREE-BAND 
SYSTEMS 



In the preceding sections, we have focused on a two-band 
tight-binding Hamiltonian, which consists of dxz and dyz or- 
bitals. However, the main idea for realizing topological insu- 
lators using two density wave order parameters with opposite 
symmetries under reflections is valid in general and applica- 
ble to more realistic multi-orbital systems. Here we extend 
our analysis to a three-band model composed of dxz, dyz, and 
dxy orbitals. In particular, we apply our idea to a more real- 
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TABLE IV: Parameters for the three-band tightbinding Hamilto- 
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istic model Hamiltonian relevant to the iron pnictide system, 
which is a representative itinerant multi-band system mani- 
festing a density wave ground state with the ordering wave 
vector Q — (tTjO).-'^'"'^^^ Here we adopt the three-orbital 
model proposed by Daghofer et al.r^ which captures the main 
physical properties of the Fe-pnictide systems. In the momen- 
tum space, the effective three-band tight-binding Hamiltonian 
is given by 



H3band - 5]^4,.(k)r'^-(k)d.,.(k), 
k,cr fi,iy 



(37) 



where 



T — 2^2 cos kx + 2ti cos ky + At^ cos kx cos ky — fic, 

T^^ — 2ti cos kx + 2^2 cos ky + At^ cos kx cos ky — Hc, 

T^^ = 2^5 (cos kx + cos ky) + Ate cos kx cos ky — fie + ^xy, 

T'^^ = T^^ =At4 sin kx sin ky , 

rpi3 ^ ^rp3iy ^ 2itTsmkx + Aits sin kx cos ky, 

y23 ^ i^rj.32y ^ 2itrsmky + Aits sin ky COS kx. (38) 



Here we use the unfolded Brillouin zone satisfying — tt < 
kx ,ky<7T as before, which corresponds to one iron atom per 
unit cell. In real iron pnictide materials, the unit cell contains 
two Fe atoms due to the buckling of the As atoms. Therefore 
the unit translations along the x (Tx) and y (Ty) directions by 
the nearest neighbor Fe-Fe distance, are not the symmetries of 
the system. However, as pointed out in Ref. 40, the system is 
invariant under the translations combined with the reflection 
Pz with respect to the xy plane, i.e. PzTx and PzTy. Then 
the eigenstates can be labeled by a pseudo-crystal momentum 
corresponding to the eigenvalues of the combined operations 
PzTx and PzTy with one iron atom per unit cell. We use this 
pseudo-crystal momentum to label states for the momentum 
space representation of the Hamiltonian in Eq. dJTJ ). 

In Eq. dJTJ l and ( l38b . //=1, 2, 3 indicate dxz, dyz, and dxy or- 
bitals, respectively. A^jj, represents the atomic potential of dxy 
orbital relative to dxz and dyz orbitals. The chemical potential 
is given by /Xc. The hopping parameters are displayed in Ta- 
ble |IVl which are determined in Ref.[33- For the 2/3 filling,— 
there are two hole pockets near the F point and two electron 
pockets at the X and Y points, which are consistent with the 
LDA calculations and ARPES measurement for LaOFeAs. 

Density wave order parameters with the momentum 
Q=(7r,0) can be described by the following Hamiltonian, 
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TABLE V: Symmetry of the charge density order parameters with 
the momentum (7r,0). Here '+' ('-') indicates 'even' ('odd') parity of 
order parameters under the corresponding symmetry operation. The 
complex off-diagonal components Dij (i 7^ j) are decomposed as 

A, = Aw + 1 Aw- 
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TABLE VI: Symmetry of the spin density order parameters with the 
momentum (7r,0). Here '+' ('-') indicates 'even' ('odd') parity of 
order parameters under the corresponding symmetry operation. The 
complex off-diagonal components Alij (i 7^ j) are decomposed as 

M,, =Mi^ + iMij. 




-ffcDW — Z^^^ {-'^Y" Dabd^^a,a'^i,b,cy, 
i.a a, 6—1 
3 

^SDW= Y. T.^-^y''^->'4,a,ai^k.AA^2- (39) 
z,cri(T2 a,fe— 1 

The order parameter represented by 3 x 3 Hermitian matrix D 
(M) has nine independent components Dij (Mij). The trans- 
formation properties of these density wave order parameters 
under the reflections PxPz and PyPz, inversion /, and time- 
reversal T are summarized in Table |V] and |Vll Notice that 
the reflections P^ and Py are not the symmetries of the sys- 
tem. The Hamiltonian is invariant only under the combined 
transformations PxPz and PyPz- 

To obtain a nodal spin density wave ground state, we con- 
sider the simplest uniform charge density wave order param- 
eter, Z)umfoim ^<iiag[do,do,dQ] with the finite diagonal com- 
ponents of Z3ii = D22 — D33 — do. This generates many 
Dirac points along an axis with a reflection symmetry in the 
momentum space, whenever a band touching occurs between 
two bands with opposite reflection parities. Now let us in- 
troduce another density wave order parameter to get a gapped 
phase. To open a full gap between neighboring bands we need 
a density wave order parameter which is odd under the reflec- 
tion symmetries PxPz and PyPz- Imposing the time reversal 
symmetry, the imaginary part of the spin density wave order 
parameter ilf/j = Im(Afi2), is the unique choice to obtain a 
topological insulator 

In Fig. [To] (a) we plot the band structure of the uniform 
charge density wave state with nonzero ^uniform 7^ along the 
ky axis. Since -Duniform preserves the PxPz reflection symme- 



FIG. 10: (Color online) The band structures of the density wave 
states with the momentum k=(7r,0) along ky axis, (a) For the uni- 
form charge density wave state with Aniform = 0.5. The locations 
of Dirac nodes are indicated by red dotted circles. The parities of 
the bands under the P.Pz symmetry are indicated by + (even) and - 
(odd), (b) A density wave phase with Duniform 7^ and Af/2 7^ at 
the same time. Here we take -Duniform=0.5 and i\//2=0.1. Each band 
is well-separated from the other bands. 
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TABLE Vll: The spin Chern numbers of the bands for a gapped 
density wave phase with Duniform 7^ and M12 7^ 0. Here we set 
Duniform = 0.5 and change the magnitude of M/2. 



try, each band has a definite reflection parity under PxPz- The 
reflection parities of the bands are also indicated in Fig.[Tol(a). 
Notice that nodal points exist between a pair of the bands with 
opposite reflection parities, which are indicated by red dotted 
circles in Fig.[l0](a). However, once we introduce a nonzero 
A//2, these nodal points disappear and a fully gapped phase 
with well-separated bands emerges. The band structure of the 
resulting gapped phase is described in Fig.[TOl(b). 

To investigate the topological property of the gapped phase, 
we compute the spin Chern numbers of the bands. Since 
the z component of the spin Sz is still conserved, the spin 
Chern number is a well-defined quantity. In Table IVIII we 
show the distribution of the spin Chern numbers for several 
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values of M(2 supporting fully gapped phases. Here Cs.n 
indicates the spin Chern number of the nth band defined as 

Cs,n = Nc^n,-\ - Nc,n,i whcTS Nc,n,-t(=-^c,n,i) denotes 
the Chern number of the nth spin-up band. We label that the 
band 1 has the highest energy and the band index n increases 
as the energy eigenvalue decreases. In the case of the gapped 
phase that is obtained by adding a small Af/2 on the nodal 
charge density wave state with -Duniform 7^ 0, only the 4th and 
5th band support nonzero spin Chern numbers shown in the 
2nd column of Table lVlH Interestingly, however, as the mag- 
nitude of M(2 increases, band gap closing and reopening oc- 
cur successively. For instance, for the uniform density wave 
state with -Dumfoim — 0.5, the first gap closing happens be- 



tween the band 1 and 2 for M/j 



0.15. As M/2 increases 



further, another fully gapped phase is obtained with the spin 
Chern numbers displayed in the 3rd column of Table lVIII It is 
interesting to notice that after the gap closing and reopening 
process, the number of the bands supporting finite spin Chern 
numbers has increased. Similar gap closing happens again for 
M/2 ~ 0.22 leading to the redistribution of the spin Chern 
numbers shown in the last column of Table I VlH Note that all 
cases with 1/3 filling give TDW insulators, while TDW phase 
with 5/6 filling occurs only for the Af/2 — 0-2 and 0.3. 



VII. SUMMARY AND DISCUSSION 

In this paper, we investigate theoretically if topological in- 
sulators can be achieved from a nodal density wave state with 
broken translational symmetry. While a nonzero density wave 
order parameter in general opens a gap between the degener- 
ate states connected by the ordering wave vector, nodal den- 
sity wave phases occur in multi-orbital systems via transla- 
tional symmetry breaking due to the distinct symmetry prop- 
erties of orbitals. Such a nodal density wave state supports a 
large number of Dirac nodes between neighboring bands. We 
have explicitly proved that a pair of inversion symmetric Dirac 
points share the same topological winding numbers in nodal 
density wave states contrary to the Dirac points in the honey- 
comb lattice. If we introduce an additional order parameter 
whose transformation property under reflection symmetries is 
opposite to that of the underlying order parameter, the system 
can be a gapped insulator at certain filling factors. Among 
those insulators, time-reversal invariant TDW insulators with 
helical edge states are identified. 

The existence of a nodal density wave ground state is 
experimentally verified in a recent ARPES measurement 



on BaFe2As2^'* and quantum oscillation experiments on 
BaFe2As2 and SrFe2As2.'*^^^ It is interesting to notice that, 
according to these experimental studies, the velocity of Dirac 
fermions is estimated to be 14 - 20 times slower than that 
in graphene.^^ This implies that the Dirac fermions in nodal 
density wave states are more susceptible to interaction effects. 
However, according to our mean field calculation, it seems 
to be difficult to realize quantum spin Hall insulators in Fe 
pnictides system, as it favors a conventional spin density wave 
state (Mq). 

Our results in general imply that transition metal materi- 
als with two-dimensional square lattice structure possessing 
partially filled t2g orbitals are good candidates for TDW in- 
sulators. In particular, in the case of the effective two-orbital 
(three-orbital) model, the 1/4 filled (1/3 filled) system is the 
most promising for the realization of TDW insulators. How- 
ever, to make a prediction on real materials with layered per- 
ovskite structure, it is important to generalize our study to 
three dimensional systems taking into account interlayer cou- 
plings. Stacking of two dimensional TDW insulators simply 
leads to a weak topological insulator'*^ Therefore identifying 
three dimensional TDW phases with a nontrivial strong topo- 
logical invariant in the layered perovskite structure is an inter- 
esting but challenging future work. 

Finally, it is worthwhile to comment the consequences of 
relaxing the constraint of time-reversal invariance. When 
an imaginary charge density wave state (D2) breaking time- 
reversal symmetry occurs in the presence of a nodal density 
wave state, a gapped topological phase with topologically pro- 
tected edge modes can be developed. In contrast to the case 
of the quantum spin Hall insulator, here the spin-up and down 
bands have the same Chern number, which gives rise to an in- 
sulator with finite Hall conductance. Interestingly, the imagi- 
nary charge density wave state is one of the competing ground 
states in iron pnictide systems, '^'^^'^^ which is expected to be 
achieved in real materials.''" Thus searching for gapped phases 
proximate to nodal density wave states is a new avenue to 
topological phases. 
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